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The phase diagram of the penetrable square-well fluid is investigated through Monte Carlo simu- 
lations of various nature. This model was proposed as the simplest possibility of combining bounded 
repulsions at short scale and short-range attractions. We prove that the model is thermodynami- 
cs ■ cally stable for sufficiently low values of the penetrability parameter, and in this case the system 
' behaves similarly to the square-well model. For larger penetration, there exists an intermediate re- 
gion where the system is metastable, with well defined fluid-fluid and fluid-solid transitions, at finite 
■ size, but eventually becomes unstable in the thermodynamic limit. We characterize the unstable 
non-extensive phase appearing at high penetrability, where the system collapses into an isolated 
blob of a few clusters of many ovelapping particles each. 
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I. INTRODUCTION 



Unlike simple fluids, complex fluids are typically characterized by a significant reduction in the number of degrees of 
freedom, in view of the hierarchy of different length and energy scales involved. As a result, coarse-grained potentials 
accounting for effective interactions between a pair of the complex fluid units adopt analytical forms that are often 
quite different from those considered paradigmatic for simple fluids^. 

An important example of this class of potentials is given by those bounded at small separations, thus indicating 
the possibility of a partial (or even total) interpenetration. This possibility, completely unphysical in the framework 
of simple fluids, becomes on the contrary very realistic in the context of complex fluids. While the true two-body 
interactions always include a hard-core part, accounting for the fact that energies close to contact raise several orders 
of magnitude, effective interactions obtained upon averaging microscopical degrees of freedom may or may not present 
this feature, depending on the considered particular system. 
^««| ' Interesting examples with no hard-core part are given by polymer solutions, where effective polymer-polymer 
interactions can be argued to be of the Gaussian former—, and star polymers and dendrimers where the so-called 
penetrable sphere (PS) model is frequently employed^— . 

In spite of their markedly different phase behaviors^, both these effective interactions have the common attributes 
of being bounded at zero separation and lacking an attractive part. The latter feature, however, appears to be 
t-H ' particularly limiting in view of the several sources of attractive interactions typical of polymer solution, such as, for 

, instance, depletion forces^, that are typically accounted through simple attractive square-well (SW) tails. 
. rH ' A tentative of combining both the penetrability at small separation and the attraction at slighty larger scale, led 
, to the introduction of the penetrable square-well (PSW) potential 8 - - — . This can be obtained either by starting from 
J-j ■ the PS model and adding an attractive well, or by starting from the SW model and reducing the infinite repulsive 
- - -' energy to a finite one. In this way, the model is characterized by two length scales (the soft core and the width of the 
well) and by two energy scales, the height e r of the repulsive barrier and the depth e a of the attractive well. 

The ratio e a / e r , hereafter simply referred to as "penetrability" , is a measure of the accessibility of the repulsive 
barrier and, as we shall see, plays a very important role in the equilibrium properties of the fluid. When e a /e r = 0, 
the PSW model reduces to the PS model (if fcsT/e r = finite, where T is the temperature) or to the SW model 
(if kBT/e a = finite). In the latter case, the model exhibits a fluid-fluid phase transition for any width of the 
attractive square well^— , this transition becoming metastable against the formation of the solid for a sufficiently 
narrow wellii. As penetrability e a /e r increases, different particles tend to interpenetrate more and more because this 
becomes energetically favorable (the precise degree depending on the e a /e r ratio). As a result, the total energy may 
grow boundlessly to negative values and the system can no longer be thcrmodynamically stable. The next question 
to be addressed is whether this instability occurs for any infinitesimally small value e a /e r > or, conversely, whether 
there exists a particular value where the transition from stable to unstable regime occurs. 

As early as the late sixties, the concept of a well-behaving thermodynamic limit was translated into a simple rule, 
known as Ruelle's criterio n 18 ' 19 , for the sufficient condition for a system to be stable. In a previous paper—, we have 
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discussed the validity of Ruelle's criterion for the one-dimensional PSW case and found that, indeed, there is a well- 
defined value of penetrability e a /e r , that depends upon the range of the attractive tail, below which the system is 
definitely stable. Within this region, the phase behavior of the fluid is very similar to that of the SW fluid counterpart. 
More recently 2 ^, we have tackled the same issue in the three-dimensional fluid. Here we build upon this work by 
presenting a detailed Monte Carlo study of the phase diagram for different values of penetrability and well width. 
In this case the PSW fluid is proven to satisfy Ruelle's criterion below a well-defined value of penetrability that is 
essentially related to the number of interacting particles for a specific range of attractive interaction. For higher 
values of penetrability, we find an intermediate region where, although the system is thcrmodynamically unstable 
(non-extensive) in the limit TV — > oo, it displays a "normal" behavior, with both fluid- fluid and fluid-solid transitions, 
for finite number of particles N. The actual limit of this intermediate region depends critically upon the considered 
temperatures, densities, and size of the system. Here the phase diagram is similar to that of the SW counterpart, 
although the details of the critical lines and point location depend upon the actual penetrability value. For even 
higher penetrability, the system becomes unstable at any studied value of N and the fluid evolves into clusters of 
overlapping particles arranged into an ordered phase at high concentration, with a phenomenology reminiscent of that 
displayed by the PS model, but with non-extensive properties. 

The remaining of the paper is organized as follows. In Section |H] we define the PSW model and in Section UTTl we 
set the conditions for Ruelle's criterion to be valid. The behavior of the system outside those conditions is studied in 
Section IIV1 where we also determine the fluid-fluid coexistence curves for the PSW model just below the threshold 
line found before; in Section [V] we determine the instability line, in the temperature-density plane, separating the 
metastable normal phase from the unstable blob phase. Section I VII is devoted to the fluid-solid transition and in 
Section [VIII we draw some conclusive remarks and perspectives. 



where e r and e a are two positive constants accounting for the repulsive and attractive parts of the potential, respec- 
tively, A is the width of the attractive square well, and a is diameter of the repulsive core. 

As discussed above, this model encompasses both the possibility of a partial intcrpenctration, with an energy cost 
typical of the soft-matter interactions given by e r , and a short-range attraction typical of both simple and complex 
fluids given by e a . Both descriptions can be clearly recovered as limiting cases of the PSW potential: For e r — > oo 
it reduces to the SW model, while for A = or e a = one recovers the PS mode l 21 ' 22 . Figure [T] displays the 
characteristics of the PSW potential (c), along with the two particular cases, SW (a) and PS (b). The interplay 
between the two energy scales e r and e a gives rise to a number of rather unusual and peculiar features that are the 
main topic of this paper. 

In order to put the PSW model in perspective, let us briefly summarize the main features of the SW and PS 
potentials. 

The SW model has a standard phase diagram typical of a simple fluid, with a fluid-fluid and a fluid-solid transitions 
in the intermediate range between the triple and the critical points in the temperature-density plane. The fluid-fluid 
transition becomes metastable, against crystallization, if the width of the well goes below a certain value that has 
been estimated to be A ps 0.25er— . 

The PS fluid, on the other hand, does not display any fluid-fluid coexistence, in view of the lack of any attractive 
interactions. The fluid-solid transition is, however, possible and highly unconventional with the formation of multiple 
occupancy crystals coupled with possible reentrant melting in the presence of a smoother repulsive interaction, such 
as a Gaussian formic. 

The PSW fluid combines features belonging to both limiting cases within a very subtle interplay between the 
repulsive and attractive energy scale that affects its thermodynamic stability*—. 



II. THE PENETRABLE SQUARE- WELL MODEL 



The PSW model is defined by the following pair potential 




(1) 



III. RUELLE'S STABILITY CRITERION 



The issue of thermodynamic stability has a long and venerable history, dating back to the late sixties^, and it is 
nicely summarized in Ruelle's textbook that is a standard reference for this problem^. 
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FIG. 1: Sketch of the PSW potential used in the present work (c). This potential interpolates between the SW potential (a) and 
of the PS potential (b). In the SW case (a), spherical particles have a perfect steric hindrance of size a (the particle diameter) 
and attractive interactions of range a + A highlighted as a halo in the picture. In the PS case (b), nearest-neighbor particles 
can partially interpenetrate, with some energy cost e r , but have no attractive tail. In the PSW there is both the possibility of 
partial interpenetration (with cost e r ) and short-range SW attraction (with energy gain e a ). 



A system is defined to be (Ruclle) thcrmodynamically stabl o 18 ' 19 if there exists a positive number B, such that for 
the total potential energy $jv for a system of N particles it holds 

<P N > —NB. (2) 

The physical rationale behind this mathematical statement is that the ratio — <3? at/TV cannot grow unboundly as N 
increases if the system is to be well behaving, but must converge to a well defined limit. This is usually referred to as 
Ruclle's stability criterion. 

Consider the PSW fluid. As density increases and temperature decreases, particles tend to lump together into 
clusters ( "blobs" ) as they pay some energetic price set by e r but they gain a (typically larger) advantage due to the 
attraction e a . Therefore, as the ratio e a /e r increases, one might expect to reach an unstable regime with very few 
clusters including a large number of significantly overlapping particles, so that $^ is no longer proportional to N. 

The ratio e a /e r ("penetrability") plays in PSW fluids a very important role, as we shall see in the following sections. 
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In Refii^ we proved that the one-dimensional (ID) PSW fluid satisfies Ruclle's criterion if e Q /e r < 1/2(1 + 1), where 
£ is the integer part of A/tr. In this case, we are then guaranteed to have a well defined equilibrium state. 

Here we show that this result can be extended to a three-dimensional (3D) PSW fluid in that Ruclle's criterion is 
satisfied if e a /^r < 1//a, where f A is the maximum number of non-overlapping particles that can be geometrically 
arranged around a given one within a distance between a and a + A. Of course, f A depends on the width of the 
attractive interaction A. For A/a < \2 — I, for instance, one has /a = 12, corresponding to a HCP closed packed 
configuration. In the following, we will use a generic d-dimensional notation and consider d = 3 at the end. 

The total potential energy of a PSW fluid formed by particles at positions ri, . . . , rjy can be written in general as 

N N 

$*(r 1 ,...,r JV ) = -£5>(|r i -r i |) (3) 

i=l jjii 

Consider now such a configuration where particles are distributed in M clusters along each direction, each made of s 
perfectly overlapped particles, and with different clusters arranged in the close-packed configuration. In the Appendix 
we prove that indeed this is the lowest possible energy configuration in the two-dimensional (2D) case. 

The total number of particles is N — M d s. As clusters are in a close-packed configuration, particles of a given cluster 
interact attractively with all the particles of those /a clusters within a distance smaller than a + A. Consequently, 
the potential energy has the form 

1 M d 
®n(M) = ^M d s (s - 1) e r - — [/a - b A (M)) s 2 e a . (4) 

The first term represents the repulsive energy between all possible pairs of particles in a given s-cluster, while the 
second term represents the attractive energy between clusters. Here 6a (M) accounts for a reduction of the actual 
number of clusters interacting attractively, due to boundary effects. This quantity clearly depends upon the chosen 
value of A/cr but we can infer the following general properties 

&a(1) = /a, 6a(M>2)</ a , lim b A (M) = 0. (5) 

M— >oo 

In the ID (with A/cr < 1) and 2D (with A/a < V% — 1) cases, $at(M) is given by Eqs. (|AT|) and (|A8|) . respectively, 
so that b A (M) = 2AJ- 1 (ID) and b A (M) = 2(4M~ 1 - M~ 2 ) (2D). In general, b A (M) must be a positive definite 
polynomial of degree d in M with no independent term, its form becoming more complicated as d increases. 
However, we need not specify the actual form of b A (M) for our argument, but only the properties given in Eq. (|5|). 
Eliminating s = N /M d in favor of M in Eq. (|4| one easily gets 

$ N (M) e r N i , . 

-e a M- d F(M) : (6) 



N 2 



where we have introduced the function 



F(M) = b A (M) - (/ A - J j . (7) 

Note that F(M) is independent of N. If e a /e r < 1/ f A , F(M) is positive definite and so §n/N has a lower bound 
(— e r /2) and the system is stable in the thermodynamic limit. Let us suppose now that e a /e r > 1//a- In that case, 
F(l) = e r /e a > but limM->oo F(M) = —(f A ~e r /e a ) < 0. Therefore, there must exist a certain finite value M = Mq 
such that F(M) < for M > M . In the ID (with A/a < 1) and 2D (with A/a < \/3 - 1) cases the values of M Q 
can be explicitly computed: 

x -1 

€ r 



Mo = [ 1 - ^ ) , (d = 1), (8) 



M, 2 ^ T 7 7 ~^ (l ^) (d = 2). (9) 

In general, it is reasonable to expect that Mo ~ (1 — e r / f A e a ) 1 . Regardless of the precise value of Mo, we have that 
lirxijv-^oo [— &n(M)]/N = oo for M > M and thus the criterion © is violated. 

This completes the proof that, if e a /e r < 1//a, the system is thermodynamically stable as it satisfies Ruelle's 
stability criterion, Eq. ©. Reciprocally, if e a /e r > 1/ f A there exists a class of blob configurations violating Eq. ©. 
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In those configurations the N particles are concentrated on a finite (i.e., independent of N) number of clusters, each 
with a number of particles proportional to N. For large N the potential energy scales with iV 2 and thus the system 
exhibits non-extensive properties. 

In three dimensions, f A = 12, 18, and 42 if A/a < \/2 - 1, y/2 - 1< A/a < \/3 - 1, and y/i - 1 < A/a < 1, 
respectively, and so the threshold values are e a /e r = 1/12, 1/18, and 1/42, respectively. There might (and do) exist 
local configurations with higher coordination numbers, but only those filling the whole space have to be considered 
in the thermodynamic limit. 

In general, Ruelle's criterion ([2"]) is a sufficient but not necessary condition for thermodynamic stability. Therefore, 
in principle, if e a /e r > 1//a the system may or may not be stable, depending on the physical state (density p and 
temperature T). However, compelling arguments discussed in RefJ^ show that the PSW system with e a /e r > 1//a 
is indeed unstable (i.e., non-extensive) in the thermodynamic limit for any p and T. Notwithstanding this, even if 
£a/e?- > V/a, the system may exhibit "normal" (i.e., extensive) properties at finite N, provided the temperature 
is sufficiently high and/or the density is sufficiently low. It is therefore interesting to investigate this regime with 
the specific goals of (i) defining the stability boundary (if any) and (ii) outlining the fate of the SW-likc fluid-fluid 
and fluid-solid lines as penetrability increases. This will be discussed in the next section, starting from the fluid-fluid 
coexistence lines. 



IV. EFFECT OF PENETRABILITY ON THE FLUID-FLUID COEXISTENCE 

We have performed an extensive analysis of the fluid-fluid phase transition of the three-dimensional PSW fluid by 
Gibbs Ensemble Monte Carlo (GEMC) simulations^—. In all cases we have started with the SW fluid (e a /e r = 0) and 
gradually increased penetrability e a /e r until disappearance of the transition. Following standard prescriptions 2 ^—, 
we construct the fluid-fluid coexistence lines using two systems (the gas and the liquid) that can exchange both volume 
and particles in such a way that the total volume V and the total number of particles N are fixed and the pressure and 
chemical potential coincide in both systems. N = 512 particles were used. By denoting with Li and Vi (i = v, I) the 
respective sizes and volumes of the vapor and liquid boxes, we used 2N particle random displacements of magnitude 
0.15Lj, N/10 random volume changes of magnitude 0.1 in ha[Vi/(V — Vi)], and N particle swaps between the gas and 
the liquid boxes, on average per cycle. 

Our code fully reproduces the results of Vega et alJ^ for the SW fluid, as further discussed below. Figure depicts 
some representative examples of the effect of penetrability on the SW results at different well widths A/ a. As A/ a 
increases, the upper limit set by Ruelle's stability condition e a /e r < 1//a decreases, and lower penetrability values 
e a /e r have to be used to ensure the existence of the transition line. In Fig. [2] values e a /e r = 1/6, 1/8, 1/11 were used 
for A/a = 0.25, 0.5, 1, respectively. Figure[5]also includes an estimate of the critical points for the PSW fluid obtained 
from the law of rectilinear diameters, as discussed in RefJ^, that is 

Pi±Pl = Pc + A (T c -T), (10) 

where pi (p v ) is the density of the liquid (vapor) phase, p c the critical density and T c the critical temperature. 
Furthermore, the temperature dependence of the density difference of the coexisting phases is fitted to the following 
scaling form 

Pi- Pg = B(T C - Tf , (11) 

where the critical exponent for the three-dimensional Ising model j3 — 0.32 was used to match the universal fluctua- 
tions. Amplitudes A and B where determined from the fit. 

A detailed collection of the results corresponding to Fig. [2ja), (b), and (c) is reported in Table HI 
Note that seemingly stable transition curves are found in all representative cases depicted in Fig. [5J thus suggesting 
a "normal" fluid behavior for the finite-size system studied. Increasing penetrability e a /e r at fixed A/ a progressively 
destabilize the transition, until a threshold value (e a /e r )th is reached where no fluid-fluid transition is observed. Upon 
changing A/a, one can then draw a line of this values in the e a /e r and A/a plane. This is depicted in Fig. [3j where 
the instability line (e a /e r )th is found to decrease as A/ a increases, thus gradually reducing the region where the fluid- 
fluid transition can be observed, as expected. The shadowed stepwise region identifies the thermodynamically stable 
region, as guaranteed by Ruelle's criterion e a /e r < 1//a discussed above. Note that points (A/a = 0.25, e a /e r = 1/6), 
(A/a = 0.5, e a /e r = 1/8), and (A/a = 1, e a /e r = 1/11), corresponding to the values used in Fig.[2]and highlighted by 
circles, lie in the 1//a < e<iAr < (eaAr)th region, that is, outside the stable range guaranteed by Ruelle's criterion. 
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FIG. 2: Fluid-fluid coexistence lines for different well widths A/a and penetrabilities e a /e r . The SW results are those by Vega 
et al.— for the same value of A/a. Circles and boxes represent the estimated critical points for the PSW and the SW fluids, 
respectively, and the dotted lines represent the coexistence curves for the PSW case, (a) A/a — 0.25 and e a /e r — 1/6; (b) 
A/a = 0.5 and e a /e r = 1/8; (c) A/a = 1 and e a /e r = 1/11. 
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A/a = 0.25, 


e a /er = 1/6 






k B T/e a 


p v cr 3 pia :i 


u v /e a 


ui/e a p v 


- fc B TlnA 3 


IH - k B T In A 3 


0.66 


0.0377(6) 0.5634(6) 


-0.343(8) 


-3.441(13) 


-2.410(7) 


-2.51(12) 


0.70 


0.0724(15) 0.5256(15) 


-0.614(16) 


-3.100(13) 


-2.253(5) 


-2.27(6) 


0.73 


0.1093(45) 0.4805(42) 


-0.862(38) 


-2.920(45) - 


-2.157(12) 


-2.29(8) 


0.75 


0.1684(95) 0.4368(95) 


-1.204(67) 


-2.682(27) 


-2.211(8) 


-2.01(2) 






A/a = 0.5, 


e a /e r = 1/8 






k B T/e a 


p v a 3 pia 3 


u v /e a 


ui/e a p v 


- fcsTln A 3 


IH - k B T In A 3 


1.00 


0.0194(4) 0.5900(7) 


-0.254(7) 


-4.687(9) 


-4.19(2) 


-4.16(5) 


1.05 


0.0319(5) 0.5841(17) 


-0.400(9) 


-4.603(14) 


-4.00(1) 


-4.01(3) 


1.10 


0.0529(8) 0.5557(8) 


-0.651(14) 


-4.365(6) 


-3.832(6) 


-3.83(4) 


1.15 


0.0799(15) 0.5173(17) 


-0.934(18) 


-4.087(15) 


-3.726(7) 


-3.76(4) 


1.20 


0.1342(37) 0.4728(40) 


-1.464(40) 


-3.777(26) 


-3.642(6) 


-3.64(2) 
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U v /t a 


ui/e a p v 


- fcsTln A 3 


III - k B T In A 3 


2.35 


0.0327(4) 0.5920(11) 


-0.693(8) 


-8.931(12) 


-8.90(2) 


-8.87(6) 


2.45 


0.0476(5) 0.5593(16) 


-1.004(11) 


-8.439(21) 


-8.66(1) 


-8.61(3) 


2.50 


0.0577(8) 0.5844(12) 


-1.201(17) 


-8.653(17) 


-8.54(2) 


-8.59(5) 


2.54 


0.0670(12) 0.5511(37) 


-1.377(25) 


-8.231(42) 


-8.48(2) 


-8.51(2) 


2.58 


0.0769(9) 0.5361(19) 


-1.556(20) 


-8.030(22) 


-8.41(1) 


-8.38(3) 



TABLE I: Vapor-liquid coexistence data from GEMC of N = 512 PSW particles with A/a = 0.25 and e a /e r = 1/6 (top table), 
A/a = 0.5 and e a /e r = 1/8 (central table) and A/a = 1.0 and e a /e r = 1/11 (bottom table). We used 10 7 MC steps. T, 
Pi, m, pi are, respectively, the temperature, the density, the internal energy per particle, and the chemical potential of the 
vapor (i = v) or liquid (i = V) phase (A being the thermal de Broglie wavelength). Numbers in parentheses correspond to the 
error on the last digits. The estimated critical points are k B T c /t a = 0.762 and p c a A = 0.307.(top table), k B T c /e a = 1.241 and 
p c a 3 = 0.307 (central table) and k B T c /e a = 2.803 and p c a 3 = 0.292 (bottom table) 



V. STABLE, UNSTABLE, AND METASTABLE PHASES 

Interestingly, in Ruelle's textbook^, the three-dimensional PSW model corresponding to point (A/a — l,e a /e r = 
1/11) is exploited as an example of "catastrophic" fluid (see especially Fig. 4 and proposition 3.2.2 both in RefJ^). 
This is clearly because this state point lies outside the stable region identified by Ruelle's criterion, as discussed. 
As already remarked, however, this criterion does not necessarily imply that outside this region the system has to 
be unstable, but only that it is "likely" to be so. There are then two possibilities. First, that in the intermediate 
region 1//a < e a Ar < (eaAr)th the system is indeed stable in the thermodynamic limit, a case that is not covered 
by Ruelle's criterion. Numerical results reported in Figs. [2] and [3] appear to support this possibility. The second 
possibility is that, even in this region, the system is strictly unstable, in the thermodynamic limit, but it appears to 
be a "normal" fluid when considered at finite N. This possibility cannot be ruled out by any simulation at finite N, 
and would be more plausible as hinted by Ruelle's arguments. 

In order to illustrate the fact that, at finite N, the system in the intermediate region 1//a < e a /er < (ea/er)th 
behaves as a normal fluid, in Fig. 0] we show two representative snapshots of the gas and the liquid phases at the point 
(A/cr = 0.5, e a /e r = 1/8) that lies just below the (e a /e r ) t h line (see Fig. [3]). In both the gas and the liquid phases, the 
structure of the fluid presents the typical features of a standard SW fluid, with no significant overlap among different 
particles. 

On the other hand, we have observed that above the threshold line (e /e r )th of Fig. [3J at a temperature close to 
the critical temperature of the corresponding SW system, the GEMC simulation evolves towards an empty box and 
a collapsed configuration in the liquid box. 

The second scenario described above can be supported or disproved by a finite-size study of the iV-dependence of 
the transition, as described below. 

Assume that at any finite the absolute minimum of the internal energy corresponds to the "collapsed" non- 
extensive configurations, referred to as "blob phase" in the following. As discussed in section Hill the internal energy 
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FIG. 3: Plot of penetrability e a /e r as a function of A/a. The displayed (e a /e r )th line separates the parameter region where 
the PSW model, with TV = 512, admits a fluid-fluid phase transition (below this line) from that where it does not. The 
shadowed stepwise line highlights the region (e a /er < 1/12 for A/a < v2— 1, £a/er < 1/18 for V2 — 1 < A/a < %/3 — 1, and 
e a /er < 1/42 for v3 — 1 < A/ a < 1) where the model is guaranteed to be thermodynamically stable for any thermodynamic 
state by Ruelle's criterion. The SW model falls on the e a /e r = axis (with finite ksT/ea). The vertical dashed arrow points 
to the SW value A/a < 0.25 below which the fluid-fluid transition becomes metastable against the freezing transition^. The 
circles are the points chosen for the calculation of the coexistence lines (Figs. [5] and O, while the crosses are the points chosen 
for the determination of the boundary phases discussed in Figs. [5] and [6] 



of these configurations scales with N 2 for large N. However, the system presumably also includes a large number 
of "normal" configurations with an internal energy that scales linearly with N . This will be referred to as "normal 
phase" . 

There is then an energy gap between the total energy associated with the normal and the collapsed configurations 
with an energy ratio of order N . For finite N and sufficiently high temperature, the Boltzmann statistical factor 
e -*iv/fc B T f co vj a p SCC i configurations (in spite of the gap) might be not sufficiently large to compensate for the 
fact that the volume in phase space corresponding to normal configurations has a much larger measure (and hence 
entropy) than that corresponding to collapsed configurations. As a consequence, the physical properties look like 
normal and one observes a normal phase. Normal configurations have a higher internal energy but also may have a 
larger entropy. If ./V is sufficiently small and/or T is high enough, normal configurations might then have a smaller free 
energy than collapsed configurations. On the other hand, the situation is reversed at larger N and finite temperature, 
where the statistical weight (i.e., the interplay between the Boltzmann factor and the measure of the phase space 
volume) of the collapsed configurations becomes comparable to (or even larger than) that of the normal configurations 
and physical properties become anomalous. This effect could be avoided only if T grows (roughly proportional to N) 
as N increases, since entropy increases more slowly with TV than <5>at. 

In a PSW fluid above the stable region (e a /e r > 1//a), we have then to discriminate whether the system is truly 
stable in the thermodynamic limit N — > oo, or it is metastable, evolving into an unstable blob phase at a given value 
of N depending on temperature and density. 

In order to shed some more light into this dual metastable/unstable scenario, we performed NVT Monte Carlo 
simulations using N = 512 particles initially distributed uniformly within the simulation box ("regular" initial con- 
dition). We carefully monitored the total potential energy of the fluid during the simulation and found that, at any 
given density, there exists a certain temperature Ti ns (p), such that the system behaves normally after 10 7 iV single 
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FIG. 4: Two GEMC simulation snapshots (TV = 512) at A/a = 0.5, e a /e r = 1/8 (below the threshold value) and k B T/e a = 1.20. 
The one on the top panel corresponds to the gas phase (p v a 3 = 0.1342), and the one on the bottom to the liquid phase 
(pia 3 = 0.4728). 



particle moves (normal phase) if T > T[ ns and collapses to a few clusters of overlapped particles (blob phase) for 

T < T ins . 

This is shown in Fig. [5] for A/a = 0.5 and two different penetrability values: e a /e r = 1/4 (upper dashed line) and 
Ea/fr = 1/7 (lower solid line). The first value lies deeply in the instability region above the threshold (e a /e r )th value 
of Fig. [31 while the second is sitting right on its top, for this value A/a = 0.5 of the well width. Also depicted are 
two snapshots of two representative configurations found under these conditions. While the particles in the normal 
phase, T > Tins; arc arranged in a disordered configuration that spans the whole box (see upper snapshot of Fig. [5]), 
one can clearly see that for T < T- ms a "blob" structure has nucleated around a certain point within the simulation 
box with a few droplets of several particles each (see lower snapshot of Fig. [5]). 

The three fluid-fluid coexistence phase diagrams displayed in Fig. [2] are then representative of a metastable normal 
phase that persists, for a given TV, up to (e a /e r )th as long as the corresponding critical point (p Cl T c ) is such that 
T c > Tins(Pe), as in the cases reported in Fig. [2] Below this instability line, the fluid becomes unstable at any density 
and a blob phase, where a few large clusters nucleate around certain points and occupy only a part of the simulation, 
is found. The number of clusters decreases (and the number of particles per cluster increases) as one moves away 
from the boundary line found in Fig. [5] towards lower temperatures. Here a cluster is defined topologically as follows. 
Two particles belong to the same cluster if there is a path connecting them, where we are allowed to move on a path 
going from one particle to another if the centers of the two particles are at a distance less than a. 

These results, while not definitive, are strongly suggestive of the fact that even the normal phase is in fact metastable 
and becomes eventually unstable in the TV — > oo limit. 

This can be further supported by a finite size scaling analysis at increasing TV, as reported in Fig. [6] in the higher 
penetrability (and hence most demanding) case e a /e r = 1/4. In obtaining these results, we used NVT simulations 
with 10 10 single particle moves in all cases. 
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FIG. 5: Regions of the phase diagram where the PSW fluid, with A/a = 0.5 and two different values of e a /e r , is expected 
to exhibit a normal phase (above the instability line) or a blob phase (below the instability line) for N = 512 particles. 
Note that the instability line corresponding to the higher penetrability case (e a /e r = 1/4, dashed line) lies above the one 
corresponding to the lower penetrability (e a /e r = 1/7, solid line). The two insets depict representative snapshots of respective 
typical configurations. 




FIG. 6: Size dependence of the instability line of Fig. [5] for the system e a /e r = 1/4 and A/a — 0.5. 
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FIG. 7: Radial distribution function for the PSW model at A/a = 0.5, ksT/t a = 1.20, and pa 3 = 0.7 for two different values 
of the penetrability parameter e a /e r '- e a /e r = 1/8 (lying below the threshold line given in Fig. [3} and e a /e r = 1/7 (that is on 
the top of it). The g(r) axis is in a log scale. 

As expected, the instability temperature line Tj ns (p) moves to higher values as N increases, at fixed density pa 3 , 
from TV = 100 to TV = 2000, and the normal phase region significantly shrinks accordingly, being expected to vanish 
in the thermodynamic limit N — > oo. 

As said before, in all the above computations we started with a regular initial condition having all particles ran- 
domly distributed in the entire available simulation box. Under these circumstances, for T < Tj ns (where all particles 
are confined into a blob of a few clusters) a large number of MC steps is required in order to find the true equilib- 
rium distribution. On the other hand, if we have a clustered configuration to start with, a much higher "melting" 
temperature Ti ns , above which one recovers a normal phase, is expected. This "hysteresis" effect is indeed observed, 
as detailed below. 

For e a /e r — 1/7, A/a = 0.5, and pa 3 — 1.0 the normal-to-blob transition occurs upon cooling at k B T/e a w 2.75. 
Inserting the obtained configuration back in the MC simulation as an initial condition, and increasing the temperature, 
we find the blob phase to persist up to much higher temperatures k B T/e a « 4. The hysteresis is also found to be 
strongly size dependent. With the same system e a /e r = 1/7, A/ a = 0.5, but for pa 3 = 0.6, we found the blob-to- 
normal melting temperatures to be k B T/e a = 2-3 for N = 256, k B T/e a = 4-5 for TV = 512, and k B T/e a = 6-7 for 
N = 1024. Analogously, in the state e a /e r = 1/4, A/a = 0.5, and pa 3 = 0.3, the results are k B T/e a = 2.1-2.2, 
k B T/e a = 3.7-3.8, k B T/e a = 9.0-9.1, and k B T/e a = 31-32 for N = 100, N = 200, N = 512, and N = 2000, 
respectively. 

In the interpretation of the size dependence of the hysteresis in the melting, one should also consider the fact that 
the blob occupies only part of the simulation box and therefore a surface term has a rather high impact on the melting 
temperature. 

Additional insights on the sudden structural change occurring on the fluid upon crossing the threshold line (e /e r )th 
can be obtained by considering the radial distribution function (RDF) g(r)3£ on two state points above and below this 
line. We consider a state point at A/a = 0.5, k B T/e a — 1.20, and pa 3 = 0.7 and evaluate the RDF at e a /e r = 1/8 
(slightly below the threshold line, see Fig. [3]) and at e a /e r = 1/7. The latter case is sitting right on the top of the 
threshold line, according to Fig. [3] The results are depicted in Fig. [7] 

Drastic changes in the structural properties of the PSW liquid are clearly noticeable. While in the normal phase 
(ea/er = 1/8) the RDF presents the typical features of a standard fluid for a soft-potential and, in particular, converges 
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to unity, in the blob phase (e a /e r = 1/7), the RDF presents a huge peak (note the log-scale) at r = and decays to 
zero after the first few peaks, a behavior that is suggestive of clustering and confinement of the system. The amplitude 
of the first maximum in the structure factor grows past the value of 2.85, which is typically reckoned for an indication 
for a freezing occurring in the system, according to Refj^ii. 

As a further characterization of the structural ordering of the system, we have also investigated a set of rotationally 
invariant local order indicators that have been often exploited to quantify order in crystalline solids, liquids, and 
colloidal gels22: 



Qi = 

where Qi m is defined as 



\ 



47T 
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m——l 



h l>2i=i N b{i)qi m {i) , , 



where N r is the number of clusters and 



1 N b (r) 

Qlm (») = ^r-y J2 Y l™ (*«) ■ ( 14 ) 
3=1 

Here Nf,(i) is the set of bonded neighbors of the i-th cluster, the unit vector fy- specifies the orientation of the bond 
between clusters i and j, and Yi m (jij) are the corresponding spherical harmonics. 

A particularly useful probe of the possible crystal structure of the system is a value of Qe close to unity (see 
Appendix A of Rcf^). Results for Qe from the PSW model are reported in Table HTl for the two values of penetrability 
considered in Fig. [5] (e a /e r = 1/4 and e a /e r = 1/7). In order to compute Qq, the center of mass of each cluster (as 
topologically defined before) is identified. Then, the cutoff distance for the nearest-neighbors "bonds" is selected to 
be approximately equal to the second minimum of g(r) (r 1.5<r). As detailed in Table [TTl we find 0.03 < Qe < 0.1 
for e a /e r = 1/4 (top table) and 0.05 < Qe < 0.12 for e a /e r = 1/7 (bottom table), depending on the considered values 
of temperature and density. These values have been computed with N = 512 particles but an increase up to N = 1024 
yields only a slight increase of Qe- Besides Qe, in Table [TTl we report other properties of the blob phases found with 
A/a = 0.5 and e a /e r = 1/4 and e a /e r = 1/7, such as the number of clusters and the internal energy per particle 
u/e a . We observe that the number of clusters is rather constant (typically 40-60) for penetrability e a /e r = 1/7. For 
the higher penetrability e a /e r = 1/4 the number of clusters is generally larger, as expected, but is quite sensitive to 
the specific density and temperature values. As for the internal energy per particle, we observe that its magnitude is 
always more than four times larger than the kinetic contribution ^ksT. 

No conclusive pattern appears from the analysis of results of Table [TTl as there seems to be no well-defined behavior 
in any of the probes as functions of temperature and density, and this irregular behavior can be also checked by an 
explicit observation of the corresponding snapshots. Nonetheless, these results give no indications of the formation of 
any regular structure. 

The final conclusion of the analysis of the fluid-fluid phase diagram region of the PSW model is that the system 
is strictly thermodynamically stable for e a /e r < 1//a and strictly thermodynamically unstable above it, as dictated 
by Ruelle's stability criterion. However, if e a /e r > 1//a there exists an intermediate region where the system looks 
stable for finite N and becomes increasingly unstable upon approaching the thermodynamic limit. 

The next question we would like to address is whether this scenario persists in the fluid-solid transition, where 
already the PS model displays novel and interesting features. This is discussed in the next section. 



VI. THE FLUID-SOLID TRANSITION 



It is instructive to contrast the expected phase diagram for the SW model with that of the PSW model. 

Consider the SW system with a width A/cr = 0.5 that is a well-studied intermediate case where both a fluid-fluid and 
a fluid-solid transition have been observed^. The corresponding schematic phase diagram is displayed in Fig. [5] (top 
panel), where the critical point is (fcsT c /e a = 1.23, p c a 3 = 0.309) in the temperature-density plane, and its triple point 
is (fcsT t /e a = 0.508, Pt<7 3 /e a — 0.00003) in the temperature-pressure plane, with pi<7 3 = 0.835 and p s a 3 = 1.28i£. 
In Refpii no solid stable phase was found for temperatures above the triple point, meaning that the melting curve 
in the pressure-temperature phase diagram is nearly vertical (see Fig. [51 top panel). Motivated by previous findings 
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TABLE II: Number of clusters, Qe parameter, and internal energy per particle for the non-extensive phases found in the case 
A/a = 0.5 and e a /e r = 1/4 (top table) and e a /e r = 1/7 (bottom table), just below the curves of Fig. [5] The parameter Qa 
was calculated on the final equilibrated particle configuration only, with a neighbor distance of 1.5cr in all cases. 



in the fluid-fluid phase diagram, we consider the PSW model with A/a = 0.5 and two different penetrability values 
£aAr = 1/15 and e a /e r = 1/8 in the intermediate region 1//a < £aAr < (ea/er)th (see Fig. [3]), where one expects a 
normal behavior for finite N, but with different details depending on the chosen penetrability. In the present case, 
the first chosen value (e a /t r = 1/15) lies very close the boundary (e a /e r = 1//a) of the stability region predicted by 
Ruellc's criterion, whereas the second chosen value lies, quite on the contrary, close to the threshold curve (e a /e r )th- 

We have studied the system by isothcrmal-isobaric (NPT) MC simulations, with a typical run consisting of 10 8 
MC steps (particle or volume moves) with an equilibration time of 10 7 steps. We considered N = 108 particles and 
adjusted the particle moves to have acceptance ratios of approximately 0.5 and volume changes to have acceptance 
ratios of approximately 0.1. Note that the typical relaxation time in the solid region is an order of magnitude higher 
than that of the liquid region. 

Consider the case e a /e r = 1/8 first. The result for the isotherm ksT/ca = 1 is reported in Fig. HI this temperature 
being smaller than the critical one fcsT c /e a = 1.241. From this figure we can clearly see the jumps in the density 
corresponding to the gas-liquid coexistence region and to the liquid-solid coexistence region. On the basis of the 
obtained results, we can foresee a phase diagram of the PSW system for this particular value of penetrability to be 
the one sketched in Fig. [5] (bottom panel). In particular, the melting curve has a positive slope in the pressure- 
temperature phase diagram, unlike the almost vertical slope of the SW counterpart, as discussed. This implies 
that penetrability allows for a "softening" of the liquid-solid transition, so the liquid and the solid can coexist at a 
temperature higher than the triple one without the need of a huge increase of pressure. 

Next we also consider a fluid with e a /e r = 1/15, just outside the Ruelle stability region, at the same temperature 
as before. The results are also reported in Fig. [8] and show no indications of a stable solid in the considered range of 
pressures, in agreement with the fact that at this very low value of penetrability the behavior of the system is very 
close to the SW counterpart. 

A specific interesting peculiarity of the PSW system in the intermediate region 1//a < £a/e?- < (eaAr)th of Fig. [3] 
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FIG. 8: Schematic phase diagram of the SW fluid for A/a = 0.5 (top panel) and phase diagram of the PSW fluid for A/a — 0.5 
and e a /e r — 1/8 (bottom panel). 



is a lack of full consistency with known thermodynamic relations. In this case, in fact, unlike the SW counterpart, 
we were unable to trace the coexistence curve between the liquid and the solid using Koike's metho d 31 ' 32 , which is 
equivalent to the numerical integration of the Clausius-Clapeyron equation 

d\nP\ Ah 1 

^r) c --jp^> p= i^t> (15) 

with Ah = hi — h s and Aw = Vi — v s , where hi and Vi denote, respectively, the molar enthalpy and volume of phase i 
(i = I for the liquid phase and i = s for the solid phase); the subscript c indicates that the derivative is taken along 
the coexistence line. Once a single point on the coexistence curve between the two phases is known one can use a 
trapezoid integration scheme^ to integrate Eq. (|T5|) . 

In our calculation, we have selected a penetrability e a /e r = 1/8 and the isotherm of Fig. HI fcsT/e a = 1, as a 
reference point. The coexistence pressure at that temperature is Pa 3 /e a w 0.475 and the molar volume jump is 
Av/a 3 ~ 1/0.775 — 1/1.313 ss 0.529. We have then calculated the molar enthalpy in the NPT ensemble by computing 
(PV + U) /N (where U is the total internal energy of the system) with the result Ah/e a -5.042 - (-7.593) = 2.551. 
Choosing a spacing in f3 of — 0.05/e a we get from Eq. (fT5"|) a predicted coexistence pressure Pa 3 /e a w 0.789 at 
ksT/ea = 1/0.95 ~ 1.053. Instead, however, at the latter temperature we found the coexistence pressure between 0.5 
and 0.6. Despite this quantitative discrepancy, Eq. (|15|l is useful to understand that the relatively mild slope of the 
PSW liquid-solid coexistence line in the pressure-temperature phase diagram is essentially due to the fact that the 
internal energies of the coexisting liquid and solid phases are not too disparate. 

A close inspection of several snapshots of the obtained solid phase suggests that, in the intermediate penetrability 
case, the obtained crystal is made of clusters of overlapping particles located at the sites of a regular crystal lattice 
with Qq sa 0.35^ and a triclinic structure characterized by a unit cell of sides a = b = c = a and angles a = j3 = 7r/3 
and 7 = cos^ 1 (l/4) (see three views of a common snapshot in Fig. [5]). 
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FIG. 9: Isotherm k B T/e a = 1 for the PSW system with A/a = 0.5 and e a /e r = 1/8 and e a /e r = 1/15, as obtained from NPT 
MC simulations with N = 108 particles. The pressure axis is in logarithmic scale. Three views of the same snapshot of the 
centers of mass of the clusters in the solid are shown on the right-hand side. 

It is worth stressing that the additional degree of penetrability, not present in the SW counterpart, is responsible 
for the coexistence of the liquid and the solid at not excessively large pressures. Clearly, we cannot rule out the 
possibility of other additional solid-solid coexistence regions at higher pressures. 

VII. CONCLUSIONS 

In this paper, we have studied the phase diagram of the three-dimensional PSW model. This model combines 
penetrability, a feature typical of effective potential in complex fluids, with a square-well attractive tail, accounting 
for typical effective attractive interactions that are ubiquitous in soft matter. It can then be regarded as the simplest 
possible model smoothly interpolating between PS (e a /e r — > 0, feT/e,. = finite) and SW (e a /e r — > 0, fcsT/e a = finite) 
fluids, as one changes penetrability e a /e r and temperature. 

We have proved that the model is thermodynamically stable when e a /e r < 1//a, as it satisfies Ruelle's stability 
criterion^. Above this value, the fluid is, strictly speaking, unstable in the thermodynamic limit, exhibiting non- 
extensive properties. For finite N, however, it displays a rather rich and interesting phenomenology. In particular, 
there exists an intermediate region 1//a < £aAr < (ca/tr )th i n the penetrability- width plane (see Fig. [3]) where the 
fluid displays normal or anomalous behavior depending on the considered temperatures and densities. For sufficiently 
large temperatures (T > Tj ns (p)) the fluid presents a metastable normal behavior with (apparently) stable liquid- 
liquid and liquid-solid transitions, provided the relative critical temperatures are above the instability line T = 7i ns . 
In this case, we have studied the effect of penetrability on the fluid-fluid transition (see Fig. [2]) close to the threshold 
line (e a /^r)th and found that in general the transition has a higher critical temperature than the SW counterpart. 
We have attributed this result to the additional degree of freedom given by penetrability that tends to oppose the 
formation of a crystal until a sufficient large density is achieved. 

Below the instability line Ti ns (p), however, different particles tend to overlap into a few isolated clusters (blobs) 
confined in a small portion of the available volume and the total energy does no longer scale linearly with the number 
of particles N. As a consequence, the fluid becomes thermodynamically unstable and its properties very anomalous 
(Fig. [SJ. The metastable region shrinks as either e a /e r or N increase (Fig. [6]). 
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FIG. 10: Sketch of a configuration with m = 4 rows and M = 6 clusters per row. 



Above the threshold line (e a /e r )th (sec Fig. [3J the fluid-fluid coexistence disappears, since in this case T; ns is too 
high to allow any phase-separation (for a given N). 

An additional interesting feature of the metastablc/unstable dualism is included in the hysteresis dependence on the 
initial condition. When the initial configuration is an unstable one (i.e., a blob) the system melts back to a normal 
phase at temperatures that are in generally significantly higher than those where the transition normal-to-blob is 
achieved upon cooling. We have attributed this behavior to the small statistical weight of the blob configuration in 
the Boltzmann sampling, in spite of its significantly larger energetic contribution. 

We have also studied the fluid-solid transition in the intermediate metastable region 1//a < e a Ar < (eaAr)th- We 
find that the solid density typically increases with respect to the corresponding SW case, due to the formation of 
clusters of overlapping particles in the crystal sites, as expected on physical grounds. The melting curve is found to 
have a relatively smooth positive slope, unlike the SW counterpart, and this anomalous behavior is also reflected in 
the thermodynamic inconsistency present in the Clausius-Clapeyron thermodynamic equation, thus confirming the 
metastable character of the phase. When penetrability is sufficiently low to be close to the Ruelle stable region, the 
system behaves as the corresponding SW system. 

One might rightfully wonder whether the finite N metastable phase presented here should have any experimental 
consequence at all. We believe the answer to be positive. Imagine to be able to craft, through a clever chemical 
synthesis process, a fluid that may be described by an effective interaction of the PSW form. Our work has then 
set the boundary for observing a very intriguing normal-to-collapsed phase by either tuning the temperature/density 
parameters, or by increasing the number of particles in the fluid. In this case, it is the finite N state, rather than the 
true thermodynamic limit N — > oo, the relevant one. 
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Appendix A: Ruelle's stability criterion in d = 2 

Let us consider the two-dimensional PSW model characterized by e a /e r and A/cr < \/3 — 1. The latter condi- 
tion implies that in a hexagonal close-packed configuration a particle can interact attractively only with its nearest 
neighbors, so that /a = 6. 

Given the number of particles N, we want to get the configuration with the minimum potential energy We 
assume that such a configuration belongs to the class of configurations described by m rows, each row made of M 
clusters, each cluster made of s perfectly overlapped particles. The centers of two adjacent clusters (in the same row 
or in adjacent rows) are separated a distance a. The total number of particles is A = rtiM s. Figure [TU] shows a sketch 
of a configuration with m — 4 rows and M = 6 clusters per row. The potential energy of an individual row is the 
same as that of the one-dimensional case^, namely 

$™w^ Ms i_! £r _ (M _ 1)s 2 ea . (Al) 
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The first term accounts for the repulsive energy between all possible pairs of particles in a given s-cluster, while the 
second term accounts for attractions that are limited to nearest neighbors if A/ a < y/3 — 1 in d = 2. The potential 
energy of the whole system is m<I> row plus the attractive energy of nearest-neighbor clusters sitting on adjacent rows 
(and taking into account the special case of boundary rows). The result is 



<&n(iti, s) = m 



Ms S —±e 



s-1 

= N——€ r - 



r - [M - 1) s 2 e a 
3m- 2 



(m - 1) [1 + 2(M - 1)] s 2 e a 



-A - (2m - l)s 



For a given number of rows m, the value of s that minimizes $Ar(m, s) is found to be 
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which is meaningful only if e a /e r > m/2(3m — 2) > 1/6. Otherwise, s*(m) = 1. Therefore, the corresponding 
minimum value is 



<f>* N (m) = 4>Ar(m, Sst (m)) 
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Let us first suppose that e a /e r < 1/6. In that case, t a je r < m/2(3m — 2) regardless of the value of m > 1 and, 
according to Eq. (|A4|) . the minimization of $^(m) is achieved with m = M = N 1 / 2 . As a consequence, Ruelle's 
stability criterion (J2J is satisfied in the thermodynamic limit with B = 3e a . 

Let us now minimize <&* N (m) with respect to m if e a /e r > m/2(3m — 2). This yields the quadratic equation 
(6 — e r /e a )m 2 — 12m + 4 = 0, whose solution is 



3 - y/3 + e r /e a 



(A5) 



The condition e a /e r > m**/2(3m** — 2) is easily seen to be equivalent to the condition e a /e r > 1/6. Therefore, the 
absolute minimum of the potential energy in that case is 
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The corresponding value of s* is 
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Comparison between Eqs. (|A5[) and (|A7j) shows that A = m^s**, i.e., the number of clusters per row equals the 
number of rows, M** = m**, as might have anticipated by symmetry arguments. 

Equation (|A6|) shows that, if e a /e r > 1/6, limAr_^ oc (— 4>^)/A = 00 and thus Ruelle's stability condition ([2]) is not 
fulfilled. 

We could have restricted to a symmetric arrangement from the very beginning, i.e., m = M and N = M 2 s, in 
which case Eq. (|A2[) yields 



$jv(M,s = N/M 2 ) 



, r 9 S — 1 

2 

N f N 
~2 



(MI 2 - 4M + 1) s 2 e a 

A 2 

M2 -lje r -(3Af 2 -4il/ + l) ^e a . 



(A8) 



The minimum value (if e a /e r > 1/6) corresponds to the value M = m»» given by Eq. (|A5|) . as expected. 
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